home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 293_02 / gradpp.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-08-23  |  11.6 KB  |  358 lines
  1. /*********************** gradpp.c ******************************
  2.  
  3.  
  4.    This program processess a set of CT scans for one
  5.    patient. It outputs two sets of files containing predefined 
  6.    views of the bone surface. One set contains distance (e.g., depth
  7.    coded surface) values and the other gradient shaded values.
  8.  
  9.  
  10.     *********************************************************/
  11.  
  12. /*************************************************************
  13.  
  14.       3-D Reconstruction of Medical Images
  15.  
  16.     Three Dimensional Reconstruction Of Medical
  17.     Images from Serial Slices - CT, MRI, Ultrasound
  18.  
  19.  
  20.    These programs process a set of slices images (scans) for one
  21.    patient. It outputs two sets of files containing nine predefined
  22.    views of bony surfaces. One set contains distance values and
  23.    the other gradient values.
  24.  
  25.    The distance values are used as 3-D spatial topographic surface
  26.    coordinate maps for geometrical analysis of the scanned object.
  27.  
  28.    The gradient values are used for rendering the surface maps on
  29.    CRT displays for subjective viewing where perception of small
  30.    surface details is important.
  31.  
  32.     Daniel Geist, B.S.
  33.     Michael W. Vannier, M.D.
  34.  
  35.     Mallinckrodt Institute of Radiology
  36.     Washington University School of Medicine
  37.     510 S. Kingshighway Blvd.
  38.     St. Louis, Mo. 63110
  39.  
  40.     These programs may be copied and used freely for non-commercial
  41.     purposes by developers with inclusion of this notice.
  42.  
  43.  
  44. ********************************************************************/
  45.  
  46. #include <stdio.h>
  47. #include <math.h>
  48. #include <ctype.h>
  49.  
  50. /* some global variables*/
  51.  
  52. int FIRSTSLICE,LASTSLICE,THRESHOLD,NLINES;
  53. float ZOOM;
  54. char DR;
  55. int  buffer[2][6][256];          /* input buffer */
  56. float fxbuf[5][256],fybuf[5][256]; /* X,Y Views floating buffers */
  57.  
  58. float huge fzbuf[256][256];        /* Z view floating buffer */
  59.  
  60. char *fnamein="ctbild.000";        /* input file name  */
  61.  
  62. succ(i)   /* (i+1) modolus 3 */
  63. int i;
  64. {
  65.     return(i==2?0:i+1);
  66. }
  67.  
  68. prev(i)   /* (i-1) modolus 3 */
  69. int i;
  70. {
  71.     return(i==0?2:i-1);
  72. }
  73.  
  74. /* Set input file - add slice number to extension */
  75. setfilename(filenum)
  76. int filenum;
  77. {
  78.     fnamein[7]=filenum/100+'0';
  79.     fnamein[8]=(filenum%100)/10+'0';
  80.     fnamein[9]='0'+filenum%10;
  81. }
  82.  
  83. /* interpolate from bottom line to top line n-1 lines  */
  84. interpolate(line,bot,top,n)
  85. int line,bot,top,n;
  86. {
  87.     int next,i,j,x,inc;
  88.     inc=top>bot?1:(-1);    /* interpolate backward or forwards ? */
  89.     next=bot+inc;
  90.     for(i=1,j=n-1;i<n;i++,j--){   /* do for each next line of interpolation */
  91.         for(x=0;x<256;x++) buffer[line][next][x]=
  92.             (buffer[line][bot][x]*j+buffer[line][top][x]*i)/n;
  93.         next+=inc;
  94.     }
  95. }
  96.  
  97. /* midpoint - fraction part of threshold transition distance */ 
  98. float midpoint(b,a)
  99. float b,a;
  100. {
  101.     return( (THRESHOLD-a) / (b-a) );
  102. }
  103.  
  104. /* get floting point distance values  */
  105. getdistances(xstart,xstop,xdir,ystart,ystop,ydir,start_slice,zdir,pass)
  106. int xstart,xstop,xdir,ystart,ystop,ydir,start_slice,zdir,pass;
  107. {
  108.     int z,x,y,i,j,start,stop,inc,line,rzoom,inter;
  109.     float remain;
  110.     FILE *fxfloat,*fyfloat,*fzfloat,*fn[2];
  111.     char filename[13];
  112.  
  113.     NLINES=0;         /* number of output lines in X,Y view directions */
  114.     remain=0;         /* remainder of interpolation after roundoff */
  115.     rzoom=ZOOM+0.5;   /* rounded zoom factor */
  116.  
  117.     sprintf(filename,"%c:xdis%d.dat",DR,pass);
  118.     fxfloat=fopen(filename,"wb");    /* X view floating output file */
  119.     sprintf(filename,"%c:ydis%d.dat",DR,pass);
  120.     fyfloat=fopen(filename,"wb");    /* Y view floating output file */
  121.  
  122.     for(i=0;i<256;i++)
  123.         for(j=0;j<256;j++)fzbuf[i][j]=256; /* clear Z view buffer */
  124.  
  125.     for(z=0;z<(LASTSLICE-FIRSTSLICE);z++){  /* For each Slice */
  126.         for(i=0;i<2;i++){               /* open next two slice files */
  127.             setfilename(start_slice+(z+i)*zdir);
  128.             fn[i]=fopen(fnamein,"rb");
  129.         } 
  130.  
  131.     inter=rzoom; /* interpolation factor assumed rounded zoom */
  132.  
  133.     /* correct interpolation factor according to floating remainder */
  134.         remain+=rzoom-ZOOM;
  135.         if(remain>=1){
  136.             inter-=1;
  137.             remain-=1;
  138.         }
  139.         else if(remain<=(-1)){
  140.             inter+=1;
  141.             remain+=1;
  142.         }
  143.  
  144.         line=0;               /* current input buffer line */
  145.  
  146.     for(j=0;j<inter;j++)  /* clear X,Y floating buffers */
  147.  
  148.         for(i=0;i<256;i++)fxbuf[j][i]=fybuf[j][i]=256;
  149.  
  150.     for(y=ystart;y!=ystop;y+=ydir){              /* For each line */
  151.  
  152.         for(i=0;i<2;i++)fseek(fn[i],(long)512*(y+1),SEEK_SET); /* skip to line*/
  153.  
  154.         fread(buffer[line][0],1,512,fn[0]);
  155.             fread(buffer[line][inter],1,512,fn[1]);
  156.             interpolate(line,0,inter,inter); /* interpolate in_between */
  157.  
  158.         for(i=0;i<inter;i++){   /* For each interplation line */
  159.                 for(x=xstart;x!=xstop;x+=xdir)  /* For each Voxel value */
  160.                     /* find threshold transition*/
  161.                     if (buffer[line][i+1][x]>=THRESHOLD){
  162.                         /* if first transition in X direction get floating distance */
  163.                         if(fxbuf[i][y]==256.0) fxbuf[i][y]=(x==xstart)?0:
  164.                         (x-xstart)*xdir-1+
  165.                             midpoint((float)buffer[line][i+1][x],
  166.                         (float)buffer[line][i+1][x-xdir]);
  167.  
  168.                         /* if first transition in Y direction get floating distance */
  169.                         if(fybuf[i][x]==256.0) fybuf[i][x]=(y==ystart)?0:
  170.                         (y-ystart)*ydir-1+
  171.                             midpoint((float)buffer[line][i+1][x],
  172.                         (float)buffer[1-line][i+1][x]);
  173.  
  174.                         /* if first transition in Z direction get floating distance */
  175.                         if(fzbuf[y][x]==256.0) fzbuf[y][x]=
  176.                             (z==0) && (buffer[line][0][x]>=THRESHOLD)?0:NLINES+i+
  177.                             midpoint((float)buffer[line][i+1][x],
  178.                         (float)buffer[line][i][x]);
  179.                     }
  180.             }
  181.             line=1-line;  /* change current input buffer line */
  182.         }
  183.         NLINES+=inter;           /* increment number of output lines */
  184.  
  185.         fclose(fn[0]);           /* close slice files */
  186.         fclose(fn[1]);
  187.         fwrite(fxbuf,1,inter*1024,fxfloat); /* write output lines to files */
  188.         fwrite(fybuf,1,inter*1024,fyfloat);
  189.         printf("did %d \n",z);
  190.     }
  191.  
  192.     fclose(fxfloat); /* close X,Y floating files */
  193.     fclose(fyfloat);
  194.  
  195.     /* write Z floating files */
  196.     sprintf(filename,"%c:zdis%d.dat",DR,pass);
  197.     fzfloat=fopen(filename,"wb");
  198.  
  199.     for(i=0;i<256;i++)fwrite(fzbuf[i],1,1024,fzfloat);
  200.     fclose(fzfloat);
  201. }
  202.  
  203. unsigned char grad(y1,y2,z1,z2,y_factor,z_factor)
  204. float y1,y2,z1,z2;
  205. int y_factor,z_factor;
  206.  
  207. {
  208.     float gx,gy,gz;
  209.     unsigned char gxint;
  210.  
  211.     /* get z and y components of gradient */
  212.     gz=(z1-z2)/y_factor;
  213.     gy=(y1-y2)/z_factor;
  214.  
  215.     /*compute gx - normalized x component of gradient */
  216.     gx=1/sqrt(1+gz*gz+gy*gy);
  217.     gxint=255*gx+0.5;      /*scale gx by 256 and round */
  218.     return(gxint);
  219. }
  220.  
  221. /* get gradient and depth shades for one image line */
  222. doline(lineg,lined,line,prevline,succline,z_factor,fg,fd)
  223.  
  224. unsigned char lineg[],lined[];  /*output buffers */
  225. int line,prevline,succline;  /* input buffer configuration */
  226. int z_factor; /* for choosing forward, backward or central differences */
  227. FILE *fg,*fd;  /* output files */
  228.  
  229. {
  230.     int i;
  231.  
  232.     /* do first pixel in line */
  233.     if(fxbuf[line][0]==256)lineg[0]=lined[0]=0;
  234.     else{
  235.         lined[0]=255-fxbuf[line][0]; /*distance shade */
  236.         lineg[0]=grad(fxbuf[line][0],fxbuf[line][1],fxbuf[prevline][0],
  237.         fxbuf[succline][0],1,z_factor);
  238.     }
  239.  
  240.     /* do rest of pixels inside line */
  241.     for(i=1;i<255;i++) if(fxbuf[line][i]==256)lineg[i]=lined[i]=0;
  242.     else{
  243.         lined[i]=255-fxbuf[line][i]; /*distance shade */
  244.         lineg[i]=grad(fxbuf[line][i-1],fxbuf[line][i+1],fxbuf[prevline][i],
  245.         fxbuf[succline][i],2,z_factor);
  246.     }
  247.  
  248.     /* do last pixel in line */
  249.     if(fxbuf[line][255]==256)lineg[255]=lined[255]=0;
  250.     else{
  251.         lined[255]=255-fxbuf[line][255]; /*distance shade */
  252.         lineg[255]=grad(fxbuf[line][255],fxbuf[line][254],fxbuf[prevline][255],
  253.         fxbuf[succline][255],1,z_factor);
  254.     }
  255.  
  256.     fwrite(lineg,1,256,fg); /* write to output files */
  257.     fwrite(lined,1,256,fd);
  258. }
  259.  
  260. /* create an gradient and distance shaded view */
  261. doviews(namedis,nameg,named,nlines)
  262.  
  263. char *namedis,*nameg,*named; /* floating file , gradiet and distance files*/
  264. int nlines;                  /* number of lines in image */
  265.  
  266. {
  267.     FILE *fg,*fd,*ffloat;
  268.     int z,i,j,k,midline;
  269.     char lined[256],lineg[256]; /* gradient and distance value buffers */
  270.  
  271.     midline=1;                  /* middle line in input buffer */
  272.     fd=fopen(named,"wb");       /* open output and input files */
  273.     fg=fopen(nameg,"wb");
  274.     ffloat=fopen(namedis,"rb");
  275.     fread(fxbuf,1,3072,ffloat); /* read first three floating lines */
  276.  
  277.     /* do first line */
  278.     doline(lineg,lined,0,0,1,1,fg,fd);
  279.  
  280.     /* do rest of lines */
  281.     for(z=0;z<(nlines-2);z++){      /*for each inside line */
  282.         doline(lineg,lined,midline,prev(midline),succ(midline),2,fg,fd);
  283.         fread(fxbuf[prev(midline)],1,1024,ffloat); /*read next floating line */
  284.         midline=succ(midline);
  285.     }
  286.  
  287.     /* do last line */
  288.     doline(lineg,lined,midline,prev(midline),midline,1,fg,fd);
  289.  
  290.     fclose(fg);  /* close all files */
  291.     fclose(fd);
  292.     fclose(ffloat);
  293. }
  294.  
  295. /**********************************************************/
  296. /**** MAIN ***** MAIN ***** MAIN ***** MAIN ***** MAIN ****/
  297. /**********************************************************/
  298. main()
  299.  
  300. {
  301.     char filename[13];
  302.     FILE *par;
  303.  
  304.     /* first get some parameters from user */
  305.     printf("Enter Zoom factor: ");
  306.     scanf("%f",&ZOOM);
  307.     printf("Enter starting scan number: ");
  308.     scanf("%d",&FIRSTSLICE);
  309.     printf("Enter ending scan number: ");
  310.     scanf("%d",&LASTSLICE);
  311.     printf("Enter Threshold value: ");
  312.     scanf("%d",&THRESHOLD);
  313.     THRESHOLD+=1024;
  314.  
  315.     DR=0;
  316.     while(!DR){
  317.         printf("Enter temporary drive: ");
  318.         DR=getchar();
  319.         if(!isalpha(DR)) DR=0;
  320.     } 
  321.  
  322.     /* get floating distance values (first pass over data)*/
  323.     printf("starting forward pass on data\n");
  324.     getdistances(0,256,1,0,256,1,FIRSTSLICE,1,1);
  325.  
  326.     /* create images */
  327.     printf("doing bottom (Z) views\n");
  328.     sprintf(filename,"%c:zdis%d.dat",DR,1);
  329.     doviews(filename,"gbo.out","dbo.out",256);
  330.     printf("doing right lateral (X) views\n");
  331.     sprintf(filename,"%c:xdis%d.dat",DR,1);
  332.     doviews(filename,"grl.out","drl.out",NLINES);
  333.     printf("doing rear (Y) views\n");
  334.     sprintf(filename,"%c:ydis%d.dat",DR,1);
  335.     doviews(filename,"gre.out","dre.out",NLINES);
  336.  
  337.     /* get floating distance values (second pass over data)*/
  338.     printf("starting backward pass on dat\n");
  339.     getdistances(255,-1,-1,255,-1,-1,LASTSLICE,-1,2);
  340.  
  341.     /* create images */
  342.     printf("doing top (Z) views\n");
  343.     sprintf(filename,"%c:zdis%d.dat",DR,2);
  344.     doviews(filename,"gto.out","dto.out",256);
  345.     printf("doing left lateral (X) views\n");
  346.     sprintf(filename,"%c:xdis%d.dat",DR,2);
  347.     doviews(filename,"gll.out","dll.out",NLINES);
  348.     printf("doing front (Y) views\n");
  349.     sprintf(filename,"%c:ydis%d.dat",DR,2);
  350.     doviews(filename,"gfr.out","dfr.out",NLINES);
  351.     printf("number of lines = %d\n",NLINES);
  352.  
  353.     /* write temporary drive and number of lines */ 
  354.     par=fopen("param.dat","wb");
  355.     fprintf(par,"%c %d",DR,NLINES);
  356.     fclose(par);
  357. }
  358.